Statistics of transmitted power in multichannel dissipative ergodic structures 
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We use the Random Matrix Theory (RMT) to study the probability distribution function and mo- 
ments of the wave power transmitted inside systems with ergodic wave motion. The results describe 
either open multichannel systems or their closed counterparts with local-in-space internal dissipa- 
tion. We concentrate on the regime of broken time-reversal invariance and employ two different 
analytical approaches: the exact supersymmetry method and a simpler technique that uses RMT 
eigenstatistics for closed non-dissipative systems as an input. The results of the supersymmetric 
method were confirmed by numerical simulation. The simpler method is found to be adequate for 
closed systems with uniform dissipation, or in the limit of a large number of weak local dampers. 



I. INTRODUCTION 

Transport through open chaotic systems is often 
viewed as a scattering process. Standard examples 
of systems of that kind are compound-nuclei, meso- 
scopic quantum dots and wires, microwave cavities 
and acoustic or ultrasonic bodies [1-12]. Incident 
waves are introduced into a disordered or irregularly 
shaped part of the structure via channels, e.g. waveg- 
uides or infinitely long ideal leads. Assuming negligi- 
ble dissipation, transport properties are obtained by 
relating incident and outgoing wave fields in terms 
of the unitary scattering matrix S. If internal dissi- 
pation is not negligible, it can be simulated by the 
action of additional open channels ||. This is the 
case, for example, with microwave cavities where non- 
perfectly reflecting walls cause loss of wave energy or 
with ultrasonic solids where internal friction acts in 
the bulk §, §, @. 

The scattering model thus applies to open systems 
both with and without internal dissipation. The scat- 
tering approach provides a useful tool for statistical 
characterization of chaotic transport. Assuming the 
wave dynamics inside the system to be ergodic so that 
the entire phase space of the system is explored [2-6], 
the scattering approach is combined with a statistical 
analysis based on Random Matrix Theory. For sys- 
tems without losses, one can make assumptions on 
the statistics of the S'-matrix J|, |H| . An alternative 
method uses a Random Matrix assumption on the 
level of the wave equation associated with the closed 
non-dissipative structure. Here the basic object is 
the Green function (resolvent) related to that wave 
equation, and the method works equally well in both 
open and closed, dissipative and non-dissipative com- 
plex structures [1-9,10,12]. Matrices from a Random 
Matrix Ensemble then replace the wave equation's 
linear differential operator, and the problem of con- 
structing various moments of the transport charac- 
teristics is expressed in terms of ensemble averages of 
the products of the resolvents. Transport character- 
istics calculated in that way are known, under certain 



conditions, to describe results of experimental mea- 
surements in systems with ergodic wave motion [1-6]. 

In this paper we are interested in characterizing the 
wave power T transmitted between a source at site i 
and a receiver at site j in a closed system with internal 
losses. The statistics of T are potentially useful for 
studies of power transmission in complex reverberant 
structures [12-14], where both mean power and the 
magnitude of its fluctuations away from the mean 
are important. In particular, we wish to calculate 
average T and T 2 with the ultimate goal to compare 
with measurements such as these of Ref. [12]. The 
complex amplitude of the transmitted wave is sim- 
ply proportional to the off-diagonal matrix element 
of the resolvent: G{E) = [E I + is I - H + iT]' 1 . 
Here H corresponds to the Hamiltonian of the closed 
non-dissipative chaotic structure. The matrix T de- 
scribes coupling to external channels or internal local- 
in-space losses, I is the identity matrix, the parameter 
e > describes uniform dissipation and E is the spec- 
tral variable. The quantity of prime interest is T = 
\dj (E)\ 2 ,i 7^ j, i.e. the product of retarded and 
advanced Greens functions (propagators): Gf^E) = 

[EI + ieI-H + iT]jl and = (Gf^E))* respec- 
tively. Except for slowly varying factors of receiver 
gain and source strength, the quantity T represents 
the ultrasonic power of Ref. [12]; see also [10,13]. 

The fluctuations in T, as measured in Ref. [12] 
were in only modest agreement with theoretical pre- 
dictions based on a simplified version of the random 
matrix approach. The moments of T were calculated 
there using a naive form of ensemble averaging, [8,12- 
14]. This relatively simple approach uses statistical 
assumptions for eigenfunctions and the real parts of 
the eigenvalues of the open (dissipative) system iden- 
tical to those of the corresponding closed system. As 
will be seen later, such an assumption is strictly justi- 
fied only for a special case of uniform dissipation. In 
a more general situation this approach fails. A proper 
treatment calls for a more elaborate technique, which 
we outline and present below. 

When losses are negligible the systems discussed in 
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Ref. [12-14] are invariant under time reversal. The 
appropriate choice for the corresponding random ma- 
trix H should therefore be a real symmetric matrix 
taken e.g. from the Gaussian Orthogonal Ensemble 
(GOE). In principle, the powerful methods of ensem- 
ble averaging we employ here can be used for such 
an ensemble, but the calculations are technically in- 
volved and will be presented in a separate publica- 
tion. 

Here we address ourselves to the somewhat simpler 
case in which H is complex Hermitian, generic for 
systems with broken time reversal invariance. Corre- 
spondingly, H is treated as a Hermitian NxN matrix 
consisting of uncorrelated centered random complex 
numbers, their variances defined by: (HijH* m ) = 
(A 2 /N) SuSjm with angular brackets indicating en- 
semble averaging. Such an if is a member of the 
Gaussian Unitary Ensemble (GUE) of random matri- 
ces [1]. Although our present results on the statistics 
of T are not directly applicable to the time-reversal 
invariant systems discussed in [12-14], they may elu- 
cidate the discrepancies found in Ref. [12] between 
measurements and the predictions of the 'naive' av- 
eraging. They also develop and illustrate the mathe- 
matical methods which will be used for a proper non- 
perturbative analysis of the time reversal-invariant 
problem. The present calculations are also relevant 
for scattering systems with broken time reversal in- 
variance as exemplified in certain chaotic billiards 
[8], optical and semiconductor superlattices [15] and 
quantum graphs [16]. In fact, our results on the dis- 
tribution of the off-diagonal elements of the resolvent 
extend earlier studies concentrated on diagonal en- 
tries for the same quantity, see [22-23] and references 
therein. Let us finally mention that there exist clear 
analogy between our research and that presented in 
the paper [19], see also the review [20]. However, 
the model considered in [19-20] did not take local 
dampers into account, but rather addressed effects of 
Anderson localization. 

The damping matrix T is in general Hermitian pos- 
itive semi-definite. In our model there is no loss of 
generality in assuming it to be diagonal. Indeed, in 
view of the rotational invariance of the Gaussian Uni- 
tary ensemble: H ^> UHU^ 1 (U -1 = W) we always 
can select the basis which diagonalizes T, bringing 
it to the form T — diag{-f,j, ...7,0, ...0}. The num- 
ber M < N of non-zero entries can be interpreted 
either as a number of equivalent open channels in the 
scattering system [3,4,9] or a number of equivalent 
localized 'dampers' in a closed system with losses. 
While we take all the 7s to be equal, the expressions 
we develop are easily generalized to the case of vary- 
ing damper strengths. It should be stressed that in 
general the matrices T and H do not commute, and 
therefore the eigenvectors and eigenvalues of the 'ef- 
fective non-Hermitian Hamiltonian' H — iT are not 
trivially related to those of H . This very fact makes 
the naive averaging incorrect. In contrast, the term 
is I interpreted as the 'uniform damping' preserves 



eigenvectors of H and just adds a uniform shift ie to 
all eigenvalues. 

The presence of an N x N random Hamiltonian H 
in the expression for the resolvent matrix G enables 
us to carry out the ensemble averaging exactly using 
the supersymmetry method [3,4,17,18,20,21]. Appli- 
cation of this non-perturbative technique leads to an 
expression for the entire probability distribution func- 
tion of T. We present the corresponding derivation 
in Section II. In Section III, we compare the results 
for the first two moments of T as obtained using the 
supersymmetry method and the methods of Ref. [12- 
14]. The results are then verified numerically by di- 
rect simulation of the model. Section IV contains 
conclusions. 



II. PROBABILITY DISTRIBUTION 

FUNCTION OF THE POWER. 
SUPERSYMMETRIC CALCULATION 

In the previous Section we defined the power T 
as a product of advanced and retarded Green func- 
tions Gij. Our goal is to compute the statistics, 
i.e. ensemble averages: (T) H , (T 2 ) h , etc., where 
subscript H designates averaging with the Gaussian 
weight exp {-^TrH^H}. At the first stage of the 
supersymmetric calculation we make use of the fol- 
lowing identities for the inverse propagator = 
[E + ie - H + iT]^ [3,4,17,18,20,21]: 

dctD^ 1 = J [dS^] [dS]exp{i£ b (E,S)}, 

detQf = (-i f J [d X *\ [rfx]exp{i£/ (£,*)}. 

Here we introduced 2 N— dimensional vectors S T = 
(SijS^) and x T = (xi jXiT); consisting of complex 
commuting or bosonic (b) variables and anticommut- 
ing or fermionic (f) variables respectively. ®b = 
diag{D,-D^} and % = diag {D, £>t} are 2Nx2N 
block diagonal matrices, and £& (E, S) = S^DbS , 
£/ (E,x) = X^/X- The negative sign of S 22 is 
necessary for convergence of the integrals in what fol- 
lows. Differentiating the first equality with respect to 
©bji an d ®fc 2 ij> an d then combining the result with 
the second equality, we obtain: 

T =£>~ 1 !D*~ 1 

= J [d$t] [d^]S* lj S li S^S 2j exp{i£(E,^)}, (1) 

where the integration involves four-component su- 
pervectors $ T = (S T ,x T ), and where £ (E,S) = 
Z b (E,S) + Z f (E, X ) = *tj>$, D = diag{D b ,® f } 
[3,4,17,18,20,21]. Because the random matrix H is in 
the exponent, T is now suitable for ensemble averag- 
ing. 

In a similar fashion, employing the Wick theorem 
one can verify the following formula necessary for the 
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calculation of an arbitrary moment of transmitted 
power (T n ) H : 

x exp (E, $)} = <S#SXS3? S%) 9 . (2) 

A shorthand notation (...)$ has been introduced 
for the 'Gaussian' integration over the supervector 
components. Hereafter we use the more convenient 
'[1,2]' ('retarded-advanced')-block notation for su- 
pervectors and supermatrices, see for example Ref. 
[4]. With the supervector * T = (Sf ,xf ,S%,xT) 
and the 4x4 supermatrices L = diag{l, 1, — 1, 1}, 
A = diag {1,1,-1,-1} the exponent in the integrand 
reads: 

£ {E, *) = EV f (I ® L) * + i* 1 (r <g> AL) * 
-* f (JL <g> L) * + ie^ f (I <g> L) 



and, Eq. (0) becomes: 



T» = (n!)- 2 (L™ ,F [¥] = S^SuS^Sij. (3) 



The ensemble averaging can now be performed 
with the aid of the identities [4]: 



(exp {itf (H ® L) = cxp |--L^rA 2 | 



AT 



4f = 4(LE^) p 4*i 4 



1/2 



StrAL = ¥ (i <8> L 1/2 AL 1/2 ) * = * f (7 ® AL) 



where we have set A = 1, and introduced a 4 x 4 
supermatrix A. Thus 



(!"% = (n!)" 2 



y [dtf] F" [*] exp jiEI't (7 <g> L) * - 4^ (T ® AL) W - ^-^A 2 - eStrAL^ , (4) 



r 



where k, I distinguish between retarded and ad- moves quartic (in ^f) term in the exponential; 2) 

vanced, (1 and 2) supermatrix blocks indices and p, ^—variables integration; 3) evaluation of the remain- 

q equal to b or / [4,17]. The next stages of super- ing integral using saddle point approximation in the 

symmetric procedure include [1,3,4,17,18,20,21]: 1) limit N — > oo. We have, after step 1): 
the Hubbard-Stratonovich transformation, that re- 

I 



(T n ) H = {niy 2 J [dR] exp j-y StrR 2 + ieNStrRA + iStrRAj 

[dtf f ] [d*] F n [*] exp {i (E&LV + (T <g> AL) *) } 



I 

Since StrRA = ^ I/'^RL 1 ' 2 ^ , for an arbitrary 4x4 supermatrix R, 
I 



(5) 



(T n ) H = {niy 2 J [dR] exp 1^-^-StrR 2 + ieNStrRA^ J [d&] [d&] F n [*] 

x exp {-i^L 1 / 2 {-Eh (g) I N - R <g> I N - iA <g> T) L 1/2 *| . (6) 

I 

Using the Gaussian nature of the integral: the following general relation, can be derived simi- 

J [d* f ] [d^]cxp{-i¥f^} = Sdctf-\ 
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larly to Eq. (|) 
[dtft] [d#] 



(7) 



xocp = (n!) a r 1 ia,» w r 1 ai> 66 J1 S , det/- 1 

Setting / = L 1 ' 2 {-Eh ®I - R®I -ik®T) L 1 ' 2 , 
we integrate out the components of \& with the help 
ofEq. g): 

(T n ) iy = /" [dR]F n [0]exp{-NC[R] +SC}, (8) 



where 



^[«]=«12,»«"21,66.« = -^4 

£ [iJ] = l^iri? 2 + Sir In (-£I 4 - i?) 
££ =ieN StrRK 
- M Str In 



J 4 - i7A {-Eh - R) 1 



See Appendix A for the details. 

{T n ) H is now written as an integral over 4x4 su- 
permatrix R. The stationarity condition for C[R], 
in the limit of large N, yields a stationary point R s , 
satisfying: R s = 1/ {—Eh — Rs)- The solution is not 
unique, it is a saddle manifold in a space of 4 x 4 su- 
permatrices, spanned by 'pseudounitary' supermatri- 
ces X: R s = - Eh/2 + im/Z^AZ = -Eh/2 - nvQ, 
where v — — E 2 / {2tt) is Wigner's semicircular 
mean density of eigenvalues (GUE, A = 1). See Refs. 
[4,17] for the explicit form of supermatrix Q. 



After integrating out local fluctuations over di- 
rections R orthogonal to the manifold of stationary 
points (the procedure is asymptotically exact for large 
N) the remaining integration goes over the manifold 
parametrized by Q: 



(T n ) H = {Tnsf n I [dQ] {Q 12 , bb Q2i,b b y 



x Sdet ~ M 
x exp{— ienuNStrQA} 



E 



h + *7^7A + inv^/QA 



(9) 



The expression for the n-th moment of power al- 
lows one to find the entire distribution function P (T) , 
cf. [17,19,20]: 

P(T) = J [dQ] 5 (T - (™) 2 Qu.hbQ21.bb) Y (Q) , 

or, for the 'scaled' power y = T/ {irv) 2 : 

V{y) = J [dQ] 8 {y - QumQ2im) Y (Q) , (10) 



where 



Y{Q) =Sdet ~ M 

x exp {— ienvNStrQA} 



E 



h + i-^lA. + inv^/QA 



Evaluation of the superintegral in Eq. ( |iTj| ) is pre- 
sented in Appendix A. The result reduces to: 



(1 i2 \ roc rl \i 



(AT~A^ eXP{ ~ e(Al " A2)} l^TATj 



(11) 

d V d V 2 J 2^/TTy (9 + V^+v) ->- 1 ^ l + y ~ X2 



d d 2 \ exp {-eVT+jj} f 1 VT+j + A 2 M 
— +y—J — ; ■m I d\ 2 7== — exp{eA 2 }(g + A 2 ) , 



I 

where 5 = (I/7 + 7) / (27w) and e = 2-kvNe. able to evaluate the remaining integral explicitly (see 

Setting the 'uniform damping' e to zero we were Appendix A for the details) and Eq. ( |ll] ) yielded: 



V{y) 



Pi 

P2 



{{9 - 1) M+1 - {9 + l) M+1 }pi + {(.9 + 1) M+1 + (g - 1) M+1 } P2 

8v /(i + y) 5 (.9 + VT+^) M+2 

= '(M + ?) } (g+(AJ + 2) ^/T+^) - (y+l)(2+(M + l)(y + 2)), 
= 2 (y + 1) (9 + [M + 2) y/l+y*) . 



(12) 
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Eqs. ([y]) and constitute the main result of this 
Section. 

At this point it is interesting to observe that Eq. 
( |l^ ) can be in fact used to cover the case of uni- 
form damping in closed system: e > 0, M = 0. 
For this we note that in the limit 7^0 (i.e. g ~ 
1/27TI/7 — > oo) and M — > oo but 7M ~ const the 
factor (g + A 2 ) / (g + X%) in the integrand of Eq. 



( |ll| ) is converted to exp {2irvjM (A2 — Ai)}. Such 
a replacement is equivalent to generating an effec- 
tive uniform damping e = 2-kv^M . The fact that the 
large number of weakly open channels (or weak local- 
in-space dampers) is essentially equivalent to uniform 
damping is well known, see e.g. [4]. Performing the 
limit g — > 00, M — > 00 when keeping the product 
2-ki/jM = e finite we find that the distribution Eq. 
( |T^ ) is reduced to: 

exp { — + y \ sinhe 



V{y) = 



(y + 1) (y + 2) - (y - 2) 1 + e^T+^j 



exp {-e^l + y} cosh e (l + eyl + y) 

+ / ; 



V( 1 + y) 3 



This distribution of transmitted power for systems 
with uniform damping is interesting and important 
on its own. 

Let us consider now a few other regimes. For the 
weakly damped system (M is fixed, j> 1), Eq. ( |l2| ) 
can be approximated by: 



v 



°\7 



(14) 



The asymptotic behavior of V (y) in the limit y 
00 for any M and g is given by: 



V{y) 



(M + 1) 

4y (M+3)/2 



{to 



nM+1 



(.9-1) 



M+ 



y(M + i)/2 



'} 

(15) 



Finally we want to compare the results of this 
Section with numerical solution of the model RMT 
problem. For this goal we numerically generate 
an ensemble of N x N Hermitian random matrices 
[H], typically choosing 1500 samples from the en- 
semble and N = 1000. The entries of the matrix 
H are constructed using a random number gener- 
ator, with (HijH* m ) — (1/N) SuSjm- To simulate 
the uniform damping and the case of a finite num- 
ber of local dampers we take T = el and T — 
diag^, 7, 7, 0, ...0} (with M < N non-zero en- 
tries) respectively. Then, for all members of our en- 
semble we generate the off-diagonal elements of the 
resolvent matrix Gij(E) = [EI + iT — H] 1 model- 
ing the response at site i due to excitation at site j, 
with E being the spectral parameter. 

We first consider the case of uniform damping: 
T = el. The modal density v for such a sys- 
tem is appro ximated by Wigner's semicircular law: 
v = \/4 — E 2 j (2ir). Therefore, for a fixed size of the 
matrices N and spectral variable E, we can explore 
the range of e by changing e. For E = the modal 
density is v — 1/ir and so we need not distinguish 
between T and y. In Fig.l, the numerically obtained 
histograms are compared with V (y) (Eq. (|ll|)) for 
several values of e. We see that numerical results 
correspond well with the theoretical curves. 




which shows that moments (y n ) exist only for n < 
(M + 1) /2. At the same time, as it follows from Eq. 
dll|), any non-zero e guarantees existence of all mo- 
ments. For large y, the asymptotic forms of the prob- 
ability distribution function at non-zero e are: 



V{y) 



\ sinh ( 



■ exp 



H 1/2 }< 



(16) 



, s o A/ cxp j-e?/ 1 / 2 ) „ 



y 



for M = 0, and for finite M respectively. 



Figure 1: Probability distribution function (Eq. 
and histograms of power for: (a) e = 1.0, (b) e = 2.0, 
(c) e = 4.0, (d) e = 6.0, as obtained numerically. Data 
are scaled to unit area. For each plot 1500 samples of 
\dj (E)\ 2 ,i / j were computed. 



This procedure was repeated for the damping ma- 
trix [r] = diaglj, 7, 7, 0, ...0} with M non-zero en- 
tries, by computing Gij(E = 0) for different combi- 
nations of parameters M and g. We The results are 
presented in Fig. 2. Again, the predictions of the su- 
persymmetry method agree well with the numerical 
results. 
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Figure 2: Probability distribution function (Eq. jl^)) and 
histograms of power for: (a) M = 2,g = 2.16, (b) M = 
6, g = 2.16 (c) M = 40, g = 400, (d) M = 40, g = 40, (e) 
M = 400, g = 400, (f ) M = 400, g = 40. Data are scaled 
to unit area. For each plot 1500 samples of \dj (E)\ 2 
were computed. We imposed the restrictions: i / j, and 
i > M, j > M for the non- uniform damping case, to avoid 
'recording' the response from damped sites or from the 
'source' site j, and to correspond to the assumptions in 
the theoretical analyses. Note, that for large values of g 
(plots (c_) and (e)), P (y) is not sensitive to either g or M 
(Eq. ©). 



III. MOMENTS OF POWER 

In this Section we analyze the first two moments of 
power y using two different approaches, both based 
on the RMT. These moments as obtained using the 
supersymmetric calculation (Eqs. (pi]), (|P2])) will be 
compared to those obtained using the naive approach. 
We first consider uniform damping and then M ^ 0. 

In the simpler, but inexact, approach Gij is con- 
structed as a modal sum, 



Gij (E) 



^ — ' E — E r — i(^ r 



and then averaged using the eigenstatistics of the un- 
damped GUE system. Here u r is the rth eigenmode 
and we call the imaginary part ( r of the eigenenergy 
E r the resonance width. 

We first consider the case when ( r is uniform: ( r = 
e, for all r [12]: 



, s2 V^Y^ U i U j U 



is E — E„ 



IE 



(18) 



On averaging Eq. (|18|) over the eigenmodes u r , and 



assuming they are uncorrelated, (y) becomes: 



(y) H 2 =E fK _ R ._ 



(E — E r — ie) (E- E r + ie) ' 

The summation over the eigenenergies E r is then re- 
placed with an integral (J2 r Nv J dE r ) 



(y) = 



-HCxD 



{Nu) dE r 



(E - E r - ie) (E-E r + ie) ' 



where u = vi — E 2 / (2ir) is the GUE modal density. 
Therefore, in the uniform damping case, the naive 
procedure produces: 

2 



(y) 



N 



2 2 



(19) 



1-KvNe 



where \ \u\ ) has been set to l/N (by normalization). 

The second moment of power is calculated in Ap- 
pendix B by means of the same approach, and is given 
by: 

(y ) = TtWNSe* + 4tt 4 i/ 4 7V% 4 

1 + 8tt 2 {Nvf e 2 - exp {-AtkNve) 

For the uniform damping case (£ r = e for all 
modes), application of the results of Section II is es- 
pecially straightforward. As already discussed, Eq. 
dill ) shows that in our model this case is realized ei- 
ther by setting M = with finite e or by letting M be 
large and 7 be small, such that 7M is finite. One can 
use Eq. ( J13| ) for this purpose, but it is more conve- 
nient to start with the first part of Eq. ( fill) . Setting 
M = and introducing: 

f = / dXt d\ 2 5 (y+l- \\) 



X 1 - X 2 . 



(A1-A2) 

we integrate by parts in Eq. (|l 

d 



exp{-e (Ai - A 2 )} , 



y n V(y)dy = -y n (nj + y? y )\F (21) 



n 

-n 2 I if 
la 



^dy = n 2 



y n - x ]dy. 



Integration with respect to y in Eq. ( |21| ) eliminates 
the delta function, and gives: 



(y) = 



/oo pi 
dXi / d\2 exp {— e (Ai 



A2)} 



Ai + A 2 
Ai — A2 



(22) 



7 



(y 2 ) 



y$dy 



(l-cxp{-2e}) 



1 



(23) 



which takes the same form as Eqs. (19) and ( |20|) 
upon substitution of 2-kvNe for e. Thus, for uniform 
damping, the 'naive' and supersymmetric methods 
agree, for both (y) and (y 2 )- This is not unexpected, 
because uniform damping with M = leaves eigen- 
statistics identical to those of closed systems, merely 
shifting all eigenenergies by is. The results ( gg ) and 
( |2g| ) can readily be reproduced by using P(y) as given 
by Eq. (p~3|) . These moments are plotted in Fig. 3 to- 
gether with the results of numerical simulations. The 
first two moments of y were obtained numerically, by 
inverting matrix EI + iT — H for each member of 
the ensemble. More precisely - we computed the col- 
umn vector Gij(E) (j is fixed, i = 1, ..N) by solving 
the algebraic equations: [SI + iT — H]Gij — Sij for a 
fixed value of E. Repeating this procedure 1500 times 
and averaging over the ensemble of H and over the 
N — 1 values of i ^ j, we obtained (y) and (y 2 ) for 
e = 1,2,4,6. As seen in Fig. 3, the correspondence 
is excellent. 
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Figure 3: (a) Log (y) and (b) Log (y 2 ) plotted as the 
functions of parameter e for the case of uniform damping. 
The solid lines represent theoretical predictions (Eqs. ( |l9| ) 
and jio])). The one sigma error bars were computed based 
on the observed variances of y and y 2 . 

Our next goal is to calculate (y) and (y 2 ) for the 
system with M equivalent dampers and without ad- 
ditional uniform damping (i.e. with e = 0). For 
this problem the probability distribution function of 
scaled power is given by Eq. (|l^ ) and closed form ex- 
pressions for the mean value of y and its variance are 



cumbersome. It is, therefore, reasonable to carry out 
the corresponding integrations numerically. In con- 
trast, the naive averaging, which now also includes 
integration over the resonance widths Cr, distributed 
[1,8] according to % 2 distribution: 



P 



M 



M 



r(M) 



M-l 



exp 



-M 



where T is mean resonance width, produces a rela- 
tively compact answer for the statistics of y (see Ap- 
pendix B): 



i \ f \2 ( v 2 2 M 



(y 2 ) 



4M 2 



e M M-l 



2 f(e M ,M) 



(24) 



(M-2) (M-l)e M 
ft e M\-l I 4Af(M-l) 

M (4Mejvf — 4ejvf — M) 



(25) 



(2M-3)e| f 



( £M + M) 4 
M 2 (2M-3). 



fl + — ) 

V M J 



4-2M 



where 6m = 271^ NT. We note that the % 2 distribu- 
tion for £ r , is strictly correct only for the case of F 
much less than mean level distance. It does, how- 
ever, correspond well with the actual distribution [4] 
for more arbitrary value of T, as long as M is large. 
In order to compare Eqs. (24) and (|25|) to corre- 



sponding results obtained by numerical integration of 
Eq. (|l2]), it is necessary to establish how T is related 
to the parameters N, M and g of the supersymmetric 
calculation. By definition, the mean scaled resonance 
width in open systems is proportional to the product 
of modal density v and average resonance width T 
[4]. Under the condition of uniform damping, the 
eigenmodes are equally damped, and the latter quan- 
tity is just equal to the individual damper strength 
e. In general, the relationship is not that simple and 
is given by Moldauer-Simonius formula [4]: 



2tt— = 2NvitT 
A 



1 M 

iVln^ 
2^ 9a 



1 



9„ 



1 

2^ 



— +la 

la 



(26) 



for the mean scaled resonance width, where A = 
1/vN is the mean eigenenergy spacing and j a is the 
coupling constant of ath channel. Note that our defi- 
nitions of 7 and eigenwidth £ r are different by factor 
of two from notation of Ref. [4]. For the uniformly 
damped system we find: 



1 



JY 



1 



1 



1 

2-KVE 



2-kveN 



(27) 



which was also shown at the end of Section II. Thus, 
we conclude that our parameter e coincides with the 
mean scaled resonance width in the limit of large 
number of equivalent weak channels. Moreover, pa- 
rameter cm in Eqs. (24), ( p5| ) is the same as 7. 

The averaging, performed in Ref. [12] for the uni- 
formly damped GOE system led to the dependence of 
the first two moments of the transmitted ultrasonic 
power on a single structural parameter, the modal 
overlap M.. M. was defined in Ref. [12] in terms of 
average imaginary part of the eigenfrequency u> r and 
modal density: M = 2n {Imuo r ) (dN/duj). We see 
that modal overlap may be identified with 7. 

(y) and (y 2 ) as predicted by supersymmetric cal- 
culation (from Ecl ([l2])) and by 'naive' averaging 
(Eqs. (24) and (|25|)) are compared in Fig. 4 and 
Fig. 5 with numerical results for several different 
values of M and 7. The prediction by the super- 
symmetry method agrees with numerical results. In 
contrast, the results of the 'naive' averaging underes- 
timate both first and second moments of the power, 
except for very large M , close to the uniform damping 
case. 

Finally, in connection with discussion of Ref. [12] 
we present the comparison of the relative variance 
(Relvar — <^y 2 ) / (y) 2 — 1) of power in Fig. 6, and 
compare the supersymmetric and naive predictions. 



IV. CONCLUSIONS 

We investigated the statistical behavior of the 
power transmitted in a closed RMT system with in- 
ternal dissipation, or an open RMT system coupled 
to the exterior via a finite number of equally strong 
channels. Using the supersymmetry method for sys- 
tems with broken time reversal invariance we derived 
an expression for the probability distribution func- 
tion for this quantity and studied its first two mo- 
ments. The theoretical predictions were compared 
to the results of numerical simulations on GUE sys- 
tems with dissipation, and to the results of a 'naive' 
theory based on the RMT eigenstatistics of a closed 
non-dissipative system. The results of the supersym- 
metric calculation agree with the numerical data for 
the full range of parameters studied. 

The naive averaging predictions are in general in- 
consistent with numerical results, because its assump- 
tions (x 2 distribution of resonance widths and decou- 
pled uncorrelated Gaussian eigenmode amplitudes) 
follow from the first order perturbation theory, valid 
for small scaled resonance width. However, because 
the x 2 distribution reduces to the exact distribution 
for the case of uniform damping or in the limit of a 
large number of weak channels, the naive theory is 
accurate in this limit, for all values of scaled level 
width. 



Acknowledgments 

This work was supported by grants from the 
National Science Foundation [CMS 99-88645 and 
CMS 0201346], by computational resources from the 
National Center of Supercomputing Applications, 
by EPSRC grant GR/R13838/01 "Random matri- 
ces close to unitary or Hermitian" and by Vice- 
Chancellor grant from Brunei University. IR would 
like to thank International Programs in Engineer- 
ing in University of Illinois at Urbana-Champaign 
for financial support, the Department of Mathemat- 
ical Sciences at Brunei University for financial sup- 
port and hospitality during his visit and Antonio M. 
Garcia-Garcia for valuable suggestions. 



Appendix A: EVALUATION OF THE 
SUPERINTEGRAL 



The result of 'Gaussian' integration in Eq. @ has 
to be simplified. To bring the integrand into a form 
convenient for saddle point integration, we use the 
series of identities for the supermatrix / = L 1 / 2 fL 1 / 2 
(L 2 = I 4 ) 

j = -Eh ® In ~ R® In ~ iA ® T 
= (l N (g>h-iT(g,(A(S- 1 ))e, (Al) 

J- 1 = (i 4 ® Jjv - (A©- 1 ) <g> iry 1 

00 

= ©- 1 ^(z) fe [(A®- 1 )0 J r] fc , 

fc=0 

where © = — EI 4 — R was introduced. Supermatrix 
/ _1 is diagonal in i and j, thus: 

f bb = L 1 / 2 f bb L 1 / 2 ,fi2 i bb ii f21,bb jj = <Sl2,bf)©21,6b- 

Substituting Sdetf^ 1 = cxp j— Str In/ j into the 

result of Gaussian integration with respect to the su- 
pervector components, and considering 

(T n ) H = Mr 2 / [dR] (^Ibb^lbbYsdetf- 1 

x exp j-y StrR 2 + ieNStr i?A j , (A2) 

we separated the terms in the exponent according to 
their order in N and obtained Eq. (||): 

C[R] = ^StrR 2 + Str\n<&, 

SC = ieNStrRK - Strln [l N -iT® (Aer 1 )] 
= ieNStr RA - MStr In [j 4 - i 7 Aer x ] . (A3) 

The last identity was proved by expanding the loga- 
rithm into the series (see Ref. [4]). 

After the Gaussian integration around the sad- 
dle point in Eq. (A2) the probability distribution 
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Figure 4: Log (y) plotted as a function of mean scaled width 7 for different number of channels: (a) M = 6; (b) M = 8; 
(c) M = 40; (d) M = 400. Here we computed G i3 (E = 0) (i / j, and i > M,j > M) for fixed i and j and averaged y 
over 1500 samples from the ensemble of H. The naive averaging prediction (solid line) is compared to the prediction 
by supersymmetry method (dotted line). The one sigma error bars were computed based on the observed variance of 

y- 



function for the scaled power is expressed as an 
integral over the manifold formed by supermarices 
Q = S^AT: 

V{y) = J [dQ]6(y-Q UM Q 21M )Y(Q), 

Q is parametrized by four commuting variables Ai, 



^2, Mi, M2 and four anticommuting a, a*, j3, (3* 
Here Ax G (l,oo), A 2 G (-1,1), and | M i| 2 = 
Af — 1, I ^2! 2 = 1 — A|. We can also introduce an- 
other set of variables, according to Ai = cosh#i , /ii = 
sinh^i expji^i}, A2 = cos#2, M2 = sin#2 exp {ifa}, 
where 9\ G (0,oo), 6 2 G (0,tt), 0i, 2 G (0, 2tt). Next 
we observe [4] 



StrQA = -2i (Ai - A 2 ) , Y (Q) = Sdet ~ M 

„2N M 



I4 + i—"/A + iirv^Qh 



' 1 + 2^X2 + Y 

1 + 27W7A1 + 7 2 



cxp{-2e7ri/iVS'tr (Ai - A 2 )} = 



ff + A; 
5 + Ai 



AI 



exp {-isirvNStrQA} (A4) 
exp {-2e-KvNStr (A a - A 2 )} , 



where 5 = (I/7 + 7) / (2-kv) , 

Qi2, b& = A*i (1 - «*«/2) (1 + W2) - a*/3/i2, 
Q 2 i,66 = M X (1 - a*a/2) (1 + /?*/?/2) + a/3> 2 , 



Q\2,bbQi\,bb = I Mi 1 2 + |Mi| 2 a*P*ap 
+ |/i 2 | 2 a*^*a/3+ |^i| 2 (^*/3-a*a) 
+ MMe^ 1+ ^a/3* (A5) 
-| M i||M2|e-^+«a*/3, 
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Figure 5: Log (t/ 2 ) plotted as a function of mean scaled width 7for different number of channels: (a) M — 6; (b) M = 8; 
(c) M = 40; (d) M = 400. (y 2 ) was obtained from the same data for Gij(E — 0) as (y) in Fig. 4. The naive averaging 
prediction (solid line) is compared to the prediction by supersymmetry method (dotted line). The error bars in cases 
(c) and (d) were computed based on the observed variance of y 2 . Error bars for cases (a) and (b) are omitted, as one 
sigma bars would misrepresent confidence intervals. Indeed, standard deviation in case (a) do not exist. 



and the integration measure was denned as 



and Grassman variables, 



dQ 



da* dp* dad/3d\i d\2 o 
(2tt(A 1 -A 2 )) 2 



Substituting Eqs. flAJ), ( |A5| ) we proceed with inte- 
gration with respect to Grassman variables. First, we 
need to expand the delta function retaining only the 
terms of zero and maximum order in these variables 
Setting Q 12 , bb Q2i ,66 = — 1 + ^i + z we 
expand the delta function, 



V(y)=S(y)- 



\ 2 



A? 



d 
dy 



dy' 



dX 2 S (y + 1 



A?) 



(Ai-A 2 )' 



exp{-e (Ai - A 2 )} 



9 + A; 
9 + Ai 



(A6) 

M 



S{y-Qf 2 Qf 1 )=S{y + l-Xl-z) = S(y + l-Xl) 
-(8>(y + l-Xl)+8>; z (y + l-Xl) (1-Xt)) 
x (A? - A2) a*f3*af3 + ... = S(y+l- Xf) 

^+»^r)*(» + l-A?)a*^ + ..., 



where we used the fact that the argument of delta 
function is linear in y (z and 1 — Af), in order to 
be able to take the differential operator out of the 
integral. Then, we calculate the integral over (f)i,<f>2 



where e = 2e-KvN and we have used 

j [dg]exp{-e(A 1 -A 2 )} 
'9 + A2 



9 + Ai 



M 



S{y + l 



A?) 



6{V), 



for the 'Efetov-Wegner' term [4,18]. Integration with 
respect to Ai is simple because of the presence of the 
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Figure 6: Relative variance of power as a function of mean scaled width 7 for different number of channels: (a) M = 6; 
(b) M = 8; (c) M = 40; (d) M = 400 The results of naive averaging (solid line) consistently underestimate the results 
by supersymmetry method (dotted line). 



delta function. 



d d 2 \ exp{-eVT+^}9(y) 
T > {y) = o{y) + [ — + y— 



d v d y 2 ) 2 % /TT^ ( g + VT+y) 

cfA 2 7= — exp {e\ 2 } (g + A 2 ) 



My) 



d_ d 2 

dy V dy 



— + y—) e (y)Z(y)- 



Using property S (y) = —yS' (y), we arrive at: 



P(y) =S(y)+6(y)y±$(y) 



(A7) 



Eq. (A.7) completes the calculation of probability 
distribution function for the scaled power. This 
equation yields Eq. ( |TT| ) upon substitution of: 
\mi y ^y-^${y) = -1. 

Finally, we set e = and derive Eq. (n2). We note 



that integral in ${y) is a table integral: 

V{v) = {dy +v W* 

x (g + i) M+1 /(g + i) -(g-i) M+1 /(g-i) 
(Af + i)(. g + Vy+T) M+1 

/(«) =2*1 (V+l,l,M + 2, 



Thus, it is possible to apply differential operator to 
obtain the final form of V (y) for this case. However, 
we notice: 

1,1 (5 + A 2 ) M 



1 



(.9 + y/y + T) } - x A2 ~ 
1 



-.d\o 



(A8) 



2\/yTT (<? + ^/yTT) 



M 



(.9 + A 2 ) M dA 2 . 



The second term in the above equation can be eval- 
uated immediately, while for the first one we can use 
an identity: 



(.9 + A 2 ) M = ]T (g + Vy + ±] 

m=0 

x^-V^+lJ m \{M-m)\ 



M—m 



Ml 
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We integrate each term in Eq. ( |A8| ) separately, 
and after the first differentiation of the result with 
respect to y, the series can be summed back, so 
that the remaining procedure becomes straightfor- 
ward and leads to Eq. (|l2|). 



Appendix B: MOMENTS CALCULATION FOR 
THE CASE OF M EQUIVALENT DAMPERS 

In this Appendix we demonstrate the intermediate 
steps leading to the Eqs. @, (§J), ©• We start 



with the modal expansion for T 2 without making an 
assumption about uniform damping: 



y 2 (tw) 



4 \ -« 



^-f E — E r — i( r E -E m + ij v 



(Bl) 



E - Ei - id E - E k + i lk 



Absence of correlation between different eigen- 
modes produces the following result for the variance 
of y: 



(y 2 ) (nv f = £ 



+E 



(E - E r - iQ r f (E-E r + i( r ) 2 



H'i \ 2 ) ( \Uj\ 2 



f£(E-E r - iQ r ) 2 (E-E l + id) 2 



{E - E r - iCr) (E - E r + t'Cr) {E - E t - i&) (E — Ei + iQ) 



(B2) 



Next, we replace summation over E r and E\ with ues in Eq. 
integration — > Nv J dE r ) and take into ac- Y 2 {ttNv (E r 
count the correlation between the GUE eigenval- 



(D2) by introducing the factor 1 



|4\ o r-KT \2 /i 



{y ) - — - — - 4 — / t —27 7 71 ,2 + 7—4 (B3J 



00 roc 



(1 - Y 2 (tzNvz)) dxdz ( N "f(\ u \ 2 ) [°° [°° (1 - Y 2 {nNvz)) dxdz 



-00 (x 2 + ( r ) ({x - zf + Cfj J-00J-00 {x-iC, r f (x- z + i(i) 2 ' 



where x =E — E r , z — E r — E\ and the Dyson two- 



level correlation function for the GUE is Y 2 (£) = 
(sin^/^) 2 . Integration over x and z in Eq. (IBS)) for 
the case of uniform damping ( r — Q — e yields Eq. 
d2Ch: 



^(l«| 4 ) 2 n 2(Nv) 2 (\u\ 2 ) 4 



(y 2 ) 



(™) 4 2C r 3 



(tti/) 



Cr + Cl r {l-Y 2 [-KNvz))dz 



2CrO 



Z 2 + {Cr+QY 



(B4) 



2 n 



+ { \u 



4 1 



, l + 8(Nv) 2 n 2 e 2 - expi-ANvTre} - ANvne , 



upon substitution (|u| 4 \/(|u| 2 \ —2, (as is the case 



for complex Gaussian random numbers) and {\u\ 
1/N. 

In the case of finite number M of weak dampers 
(e = 0) the ensemble averaging includes an integra- 
tion with over a distribution of widths, given by [1,8]: 

M-i 



:ca m m fc, 



- . 1 — 1 exp <j -Adkr \ , (B5) 

r ) r (m) V r j \ r J 

where T is average resonance width and T (M) is a 
Gamma function. Starting with Eq. ( pl| ) we average 
y over the eigenmodes and eigenenergies to get: 



(y) {n V ) 2 = Nv(\u\ 2 ) 2 J 



(B6) 



Which becomes Eq. upon integration overp (£ r ) 
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We indicate this averaging by overbar: 

, 2 ttNv M 



r m - 1 



To evaluate the second moment of power we integrate 
over x, and then over ( r and Q in Eq. (B3): 



r Cr+0 f° (l-y 2 (7TiV^))^ 



2CrO J-oo Z 2 + (C r + 0)^ 



(B7) 



The remaining integral with respect to z in J it 
is convenient to use the Fourier transform of Y% (£), 
which has a form: 

/oo 
y 2 (Oexp{27ri£g}d£ = l-M,M < 1, 
-OO 

6(g)=0,|?|>l, 



(B8) 



Cr + C 
' 4C 



~t f_ x i 1 H) cxp { ~ Nv {Cr + Ci)} dq - 



The average over ( r and Q is now straightforward. 
Finally, substituting the result into the Eq.(|B7|), and 



taking into account / /\M / = 2, we obtain 



the second moment in its closed form (Eq. (|2J 
The derivation presented above assumes that reso- 
nance widths Cr and eigenmodes u r are statistically 
independent. 
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